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ABSTRACT 

In this era of precision cosmology, a detailed physical understanding on the 
evolution of cosmic baryons is required. Cosmic magnetic fields, though still 
poorly understood, may represent an important component in the global cosmic 
energy flow that affects the baryon dynamics. We have developed an Eulerian- 
based cosmological magnetohydrodynamics code (CosmoMHD) with modern 
shock capturing schemes to study the formation and evolution of cosmic struc- 
tures in the presence of magnetic fields. The code solves the ideal MHD equations 
as well as the non-equilibrium rate equations for multiple species, the Vlasov 
equation for dynamics of collisionless particles, the Poisson's equation for the 
gravitational potential field and the equation for the evolution of the intergalac- 
tic ionizing radiation field. In addition, a detailed star formation prescription 
and feedback processes are implemented. Several methods for solving the MHD 
by high-resolution schemes with finite-volume and finite-difference methods are 
implemented. The divergence-free condition of the magnetic fields is preserved 
at a level of computer roundoff error via the constraint transport method. We 
have also implemented a high-resolution method via dual-equation formulations 
to track the thermal energy accurately in very high Mach number or high Alfven- 
Mach number regions. Several numerical tests have demonstrated the efficacy of 
the proposed schemes. 

Subject headings: cosmology: theory — magnetohydrodynamics — methods: 
numerical — shock waves 



1. Introduction 

While the concordance cold dark matter cosmological model (ACDM; Krauss & Turner 
1995; Ostriker & Steinhardt 1995; Bahcall et al. 1999; Spergel et al. 2006) has been shown 
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to be remarkably successful in many respects, there remain several apparent discrepancies 
between the model and the real Universe. For example, the centers of simulated X-ray 
clusters of galaxies do not agree with observations, manifested in their inability to match even 
the simplest relations such as the temperature-luminosity relation (e.g., Ponman, Cannon, & 
Navarro 1999), among others. Closely related to this problem is the so-called "cooling flow" 
problem in the cores of many clusters, where the cooling time is less than the Hubble time 
(e.g., Fabian 1994; White, Jones, & Forman 1997; Peres et al. 1998; Allen 2000) and, hence, 
in the absence of heating sources, the gas will cool and flow towards the center. This is 
confirmed by simple hydrostatic models (Fabian & Nulsen 1977; Cowie & Binney 1977) and 
simulations (Suginohara & Ostriker 1998; Yoshida et al. 2002), with conventional physics, 
i.e., gravity and cooling. What is most puzzling is that, observationally, this supposedly 
cooling X-ray gas shows up neither in soft X-rays nor in other cooler forms in the expected 
amounts (e.g., Peterson et al. 2003; Kaastra et al. 2004). Solutions such as supernova 
heating appear incapable of rectifying the problem (Binney 2000). 

Based on a polytropic model taking into account both a removal of low entropy gas 
due to star formation and an addition of feedback energy, Ostriker et al. (2005) show that 
if the feedback energy is of the order of a few percent of the rest mass of Supermassive 
black holes, a consistent X-ray cluster model can be constructed. Supermassive black holes 
(SMBH) are now strongly believed to lurk at the centers of probably every massive galaxy 
(Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 
2002). One of the most prominent known forms of feedback from the growth of SMBH that 
is mechanically tightly coupled to surrounding gas is the powerful radio jets and giant radio 
lobes, whose energy is often observed to exceed 10^~^ MQC^ergs (Urry 2000; Kronberg et al. 
2001). Therefore, there is a good reason to believe that the inclusion of Active Galactic Nuclei 
(AGN) feedback, in form of radio jets/lobes, in direct simulations may bring theoretical 
predictions into good agreement with observations. In this case, our understanding of X-ray 
cluster formation will be much more complete and physically sound, which in turn would 
allow us to remove systematic uncertainties with regard to determinations of dark matter and 
dark energy, among others, afforded by accumulating X-ray data and upcoming Sunyaev- 
Zel'dovich (SZ) cluster data. Clearly, it has become very pressing to properly model AGN 
feedback in cosmological simulations. 

In this era of precision cosmology, multiple, independent constraints on cosmological 
parameters are vital. Weak gravitational lensing of galaxies by large-scale structure provides 
one of the most powerful probes of matter distribution in the Universe (e.g., Refregier 2003). 
Since weak lensing directly probes the mass distribution, which is dominated by dark matter, 
it is often assumed that baryonic physics is not that important. With powerful future weak 
lensing surveys such as PanSTARRS, SNAP and LSST, which have the potential capability 
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of reducing statistical errors to about 1% level with regard to power spectrum and dark 
energy equation-of-state determinations, it will be necessary to understand hitherto "less 
important" effects such as the physics of cosmic baryons. Baryons constitute a significant 
fraction (~ 17%) of total mass in the Universe (Spergel et al. 2006) and could make a 
significant contribution to the total matter power spectrum. On sub-megaparscc scales, 
baryons are known to be subject to different physics than the dark matter, and important 
differences between distributions of baryons and dark matter on small scales (< lOOkpc) 
are well observed, such as in our own galaxy. It is, however, often assumed that baryons 
follow dark matter on large scales (> IMpc), where signals of weak gravitational lensing by 
large-scale structure are detected to in turn determine the total matter distribution on these 
large scales. Estimates based on observed radio-lobe energy from AGN suggest that a large 
volume fraction of the intergalactic medium may be moved and filled with magnetic bubbles 
(Furlanetto & Locb 2001; Kronbcrg et al. 2001; Levine & Gnedin 2005). As a result, the 
effect of large-scale movement of baryonic gas under the infiuence of giant radio lobes/jets 
from AGN on the total matter power spectrum might be as large as 30% based on a simplified 
spherical model (Levine & Gnedin 2006). While the exact effect on the intergalactic medium 
of giant radio jet/lobes from AGN is presently highly uncertain due to lack of any adequate 
treatment, these initial rough estimates of potential effects point to the need of much more 
detailed investigations. It is prudent to re-examine with greater care the assumption that 
baryons follow dark matter on large scales in light of the potential capabilities of future 
weak lensing surveys. Cosmological parameters determined with claimed achievable high 
statistical accuracy based on gravitational weak lensing can be realized only after systematic 
effects on baryons due to powerful AGN feedback can be convincingly demonstrated and 
removed. Magnetohydrodynamics simulations provide the necessary tools to tackle this 
important problem. 

The idea of AGN feedback is not new (e.g., Silk & Rccs 1998; Kronbcrg et al. 2001; 
Churazov et al. 2001; Quihs, Bower, & Balogh 2001; Bruggen & Kaiser 2001; Binney 2004; 
Omma & Binney 2004; Ruszkowski, Bruggen & Begelman 2004; Henrik et al. 2004; Scan- 
napieco & Oh 2004; Begelman 2004; Dalla Vecchia et al. 2004; Begelman & Ruszkowski 
2005; Bruggen, Ruszkowski, & Hallman 2005; Croton et al. 2006). Recent SPH (smoothed- 
particle-hydrodynamics) simulations have taken a major step to implement such an AGN 
feedback (Hopkins et al. 2005a,b,c,d,2006; Robertson et al. 2006; Sajacki & Springel 2006) 
through a deposition of thermal energy in the central regions of galaxies where massive gas 
accretion onto the central SMBH takes place. The effects are clearly demonstrated to be 
significant, often with dramatic effects, such as a complete sweep-out of interstellar medium. 
But real MHD modeling is largely lacking. Undoubtedly radio jets are highly anisotropic 
and the overall dynamics and thermodynamics of magnetic jets/bubbles and baryons on 
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scales ranging from cluster cores to large-scale structure requires solving the propagation of 
relativistic jets using cosmological MHD codes in a complex, dynamic, cosmological setting. 
As an example, a highly collimated relativistic radio jet may not deposit most of its energy 
in the immediate neighborhood, whereas a thermal and spherical deposition of the energy, 
as has been implemented in SPH simulation, would cause an interaction with the immediate 
surrounding gas at the outset. 

It has been well established that there are widespread magnetic fields in the intracluster 
medium (see Carilli & Taylor 2002; Govoni 2006 for recent reviews). There has also been 
tantalizing evidence for magnetic fields in the wider intergalactic medium (IGM) (Kim et 
al. 1989). The origin of these large-scale magnetic fields is still unknown though several 
suggestions have been made (e.g., Kulsrud et al. 1997; Furlanetto & Loeb 2001; Kronberg et 
al. 2001). Furthermore, it is presently mostly unknown whether these magnetic fields have 
played an important role or not, in systems ranging from large-scale structure formation, 
to galaxy clusters, to cluster core regions and to galaxy formation itself. There have been 
some studies aimed at investigating the role of magnetic fields in the intracluster (ICM) 
and IGM (Ryu et al. 1998; Dolag et al. 2002; Dolag et al. 2005; Briiggen et al. 2005). 
Numerically, substantial efforts have gone into developing MHD solvers with SPH (Dolag 
et al. 1999, 2002), though the divergence-free condition for magnetic fields is still difficult 
to handle (Ziegler et al. 2006). For a grid-based Eulerian approach, studies by Ryu et al. 
(1998) and Briiggen et al. (2005) have included the magnetic fields passively, i.e., they have 
omitted the feedback of magnetic fields to the medium (both in the momentum equation 
and energy equation). Consequently, we are not aware of any Eulerian-based cosmological 
simulations where magnetic fields are included self-consistently. 

In this paper, we present a highly accurate cosmological MHD code, named CosmoMHD, 
for cosmological simulations that involve magnetic fields. The outline of our paper is as 
follows. We first present an overview of the CosmoMHD code and various physics packages 
included in CosmoMHD in ^ We then describe the basic solvers for ideal MHD adopted in 
our code in ^ Numerical tests are presented in ^ Further discussions are given in ^ 

2. CosmoMHD Code and Physics Packages 

2.1. CosmoMHD Code 

The framework of our CosmoMHD code is based on the combination of the TVD-ES 
cosmological hydrodynamics (HD) code (Ryu et al. 1993) and an MHD code on an adaptive 
mesh refinement (AMR-MHD) grid (Li & Li 2003). We have integrated the cosmological 
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solvers for dark matter, atomic physics and star formation feedback processes in the TVD- 
ES code with the MHD solvers in AMR-MHD code for gas dynamics (AMR features are 
not utilized in CosmoMHD). The hydro TVD-ES code has been extensively used for cosmo- 
logical applications, including Lyman-alpha forest (Cen et al. 1994), warm- hot intergalactic 
medium (Cen & Ostriker 1999a), damped Lyman-alpha systems (Cen et al. 2003), cosmo- 
logical chemical evolution (Cen & Ostriker 1999b; Cen, Nagamine, & Ostriker 2005) and 
galaxy formation (Nagamine et al. 2000,2001a,b,2004,2005a,b,2006). The MHD code has 
been extensively used for astrophysical jet simulations (e.g., Li et al. 2006, Nakamura, Li, 
& Li 2006). 

Consistent with the notations in Ryu et al. (1993), the basic equations for CosmoMHD 
in comoving coordinates are 

^ + -V-(pv) = (1) 
at a 

^^^ + -V-(pvv+p*-BB) = --pv--pV<l' (2) 
ot a a a 

^ + -V-[(E + p>-(B.v)B] = -^i?- V-V$ + ^^, (3) 

^-lvx(vxB) = -Ab (4) 
ot a 2a 

where we have implicitly assumed 7 = 5/3. Here, p is the comoving density, v is the proper 

peculiar velocity, p is the comoving gas pressure, B'^/2 is the comoving magnetic pressure, 

p* = p + B'^/2 is the total comoving pressure, E = ^pv^ + |p + ^B"^ is the total peculiar 

energy per unit comoving volume, $ is the proper peculiar gravitational potential from both 

dark matter and self-gravity, t is the cosmic time, and a is the expansion parameter. Note 

that, from the induction equation above [Eq. (jlj)], the magnetic field in the comoving frame 

expands as -Bcomoving oc a~^/^. Since -Bcomoving = -Bproperfl^''^, the magnetic field in the proper 

frame expands as -Bpropcr a~^, consistent with the flux conservation. 

In addition to the total energy equation, we have implemented solvers for two auxiliary 
equations: the modified entropy equation given in Ryu et al. (1993) and the internal energy 
equation in Bryan et al. (1995). This implementation allows us to perform better compar- 
isons among different approaches in tracking the thermal energy accurately. The equations 
are: 

95 1 ,^ , 2d ^ 

+ _V.(5v) = S, (5) 

at a a 

UPC 1 2(2 T) 

^ + -V-(pev) = pe + ^^V-v, (6) 

Ot a a a 

where S = p/ p"'^^ is the comoving modified entropy, and e is the internal energy. Note that 
Eqs. ([5]) and ([6]) remain the same even when the magnetic fields are present. 



- 6 - 



The CosmoMHD code integrates five sets of equations simultaneously: the ideal MHD 
equations for gas dynamics and magnetic field dynamics [Equations ([1]) to Qj], rate equa- 
tions for multiple species of different ionizational states (including hydrogen, helium and 
oxygen), the Vlasov equation for dynamics of coUisionless particles, the Poisson's equation 
for obtaining the gravitational potential field and the equation governing the evolution of the 
intergalactic ionizing radiation field, all in cosmological comoving coordinates. The MHD 
code consists of several approaches for solving the MHD (or HD) by high-resolution schemes 
with finite-volume and finite-difference methods. The code preserves conservative quanti- 
ties and the divergence-free condition of the magnetic fields. The code is fully parallelized 
with OPENMP directives. The MHD solver can also be used as an hydro solver when the 
magnetic fields are zero. We will describe the MHD solvers in more details in ^ The 
rate equations are treated using sub-cycles within a hydrodynamic timestep due to much 
shorter ionization timescales. Dark matter particles are advanced in time using the standard 
particle-mesh (PM) scheme. The gravitational potential on a uniform mesh is solved using 
the Fast Fourier Transform (FFT) method. 

The simulation computes for each cell and each timestep the detailed cooling and heating 
processes due to all the principal line and continuum processes for a plasma of primordial 
composition (Cen 1992). Metals ejected from star formation are followed in detail in a 
time-dependent, non-equilibrium fashion. Cooling due to metals is computed using a code 
based on the Raymond-Smith code assuming ionization equilibrium (Cen et al. 1995): at 
each timestep, given the ionizing background radiation field, we compute a lookup table 
for metal cooling in the temperature-density plane for a gas with solar metallicity, then 
the metal cooling rate for each gas cell is computed using the appropriate entry in that 
plane multiplied by its metallicity (in solar units). In addition, Compton cooling due to the 
microwave background radiation field and Compton cooling/heating due to the X-ray and 
high energy background are also included. 

We follow star formation using a well defined (heuristic but plausible) prescription used 
by us in our earlier work (Cen & Ostriker 1992,1993) and similar to that of other investigators 
(Katz, Hernquist, & Weinberg 1992; Katz, Weinberg, & Hernquist 1996; Steinmetz 1996; 
Gnedin & Ostriker 1997). A stellar particle of mass m* = c^^m^^At/t^, is created (the same 
amount is removed from the gas mass in the cell), if the gas in a cell at any time meets the 
following three conditions simultaneously: (i) flow contracting, (ii) cooling time less than 
dynamic time, and (iii) Jeans unstable, where At is the timestep, = max(tdyn, 10''yrs), 
tdyn is the dynamical time of the cell, mgas is the baryonic gas mass in the cell and c* ~ 0.07 
is star formation efficiency. In essence, we follow the classic work of Eggen, Lynden-Bell & 
Sandage (1962) and assume that the dynamical free-fall and galaxy formation timescales are 
simply related. Each stellar particle has a number of other attributes at birth, including 
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formation time ti, initial gas metallicity and the free-fall time in the birth cell tdyn- The 
typical mass of a stellar particle in the simulation is about one million solar masses; in other 
words, these stellar particles are like globular clusters. 

Stellar particles are subsequently treated dynamically as coUisionless particles. But 
feedback from star formation is allowed in three forms: ionizing UV photons, supernova 
kinetic energy in the form of galactic siiperwinds (GSW), and metal-enriched gas, all being 
proportional to the local star formation rate. The temporal release of all three feedback com- 
ponents at time t has the same form: f{t,ti,tdyn) = {^/tdyn)[{t-ti)/tdyn] exp[-(t - ti)/trfj/„]. 
Within a timestep rft, the released GSW energy to the IGM, ejected mass from stars into 
the IGM and escape UV radiation energy are ecsw f{t)ti-,idyn)'fn*c^dt, e„iassf {t,ti-,tdyn)i^*dt 
and fesc{Z)euv{Z)f{t, U, tdyn)Tn^c^dt. We use the Bruzual-Charlot population synthesis code 
(Bruzual & Chariot 1993; Bruzual 2000) to compute the intrinsic metaUicity-dependent UV 
spectra from stars with Salpeter IMF (with a lower and upper mass cutoff of 0.1 Mq and 
125 Mq). Note that euv is no longer just a simple coefficient but a function of metallicity. The 
Bruzual-Charlot code gives euv = (1.2 x 10"^ 9.7x 10-^ 8.2 x IQ-^ 7.0 x 10"^ 5.6x 10-^ 3.9 x 
10-^ 1.6 X 10"^) at Z/ Zq = (5.0 x 10-^ 2.0 x 10"^, 2.0 x 10-\ 4.0 x 10-\ 1.0, 2.5, 5.0). We 
also implement a gas metallicity-dependent ionizing photon escape fraction from galaxies 
in the sense that higher metallicity hence higher dust content galaxies are assumed to al- 
low a lower escape fraction; we adopt the escape fractions of fesc = 2% and 5% (Hurwitz 
et al. 1997; Deharveng et al. 2001; Heckman et al. 2001) for solar and one tenth of solar 
metallicity, respectively, and interpolate/extrapolate using a linear log form of metallicity. 
In addition, we include the emission from quasars using the spectral form observationally 
derived by Sazonov, Ostriker, & Sunyaev (2004), with a radiative efficiency in terms of stel- 
lar mass of cqso = 2.5 x 10^^ for hu > 13.6eV. Finally, hot, shocked regions (like clusters 
of galaxies) which emit ionizing photons due to bremsstrahlung radiation are also included. 
The UV component is simply averaged over the box, since the light propagation time across 
our box is small compared to the timesteps. The radiation field (from leV to lOOkeV) is 
followed in detail with allowance for self-consistently produced radiation sources and sinks 
in the simulation box and for cosmological effects, i.e., radiation transfer for the mean field 
Ji, is computed with stellar, quasar and bremsstrahlung sources and sinks due to Lya clouds 
etc. In addition, a local optical depth approximation is adopted to crudely mimic the lo- 
cal shielding effects: each cubic cell is flagged with six hydrogen "optical depths" on the 
six faces, each equal to the product of neutral hydrogen density, hydrogen ionization cross 
section and scale height, and the appropriate mean from the six values is then calculated; 
equivalent ones for neutral helium and singly-ionized helium are also computed. In comput- 
ing the global sink terms for the radiation field, the contribution of each cell is subject to the 
shielding due to its own "optical depth" . In addition, in computing the local ionization and 
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cooling/heating balance for each cell, we take into account the same shielding to attenuate 
the external ionizing radiation field. 

Galactic superwinds energy and ejected metals are distributed into 27 local gas cells 
centered at the stellar particle in question, weighted by the specific volume of each cell. 
We fix Cmass = 0.25. GSW energy injected into the IGM is included with an adjustable 
efficiency (in terms of rest-mass energy of total formed stars) of ecsw, which is normalized 
to observations for our fiducial simulation with ecsw = 3 x 10~^. If the ejected mass 
and associated energy propagate into a vacuum, the resulting velocity of the ejecta would 
be {easw/^massY^'^c = 1469km/s. After the ejecta has accumulated an amount of mass 
comparable to its initial mass, the velocity may slow down to a few hundred km/s. We 
assume this velocity would roughly correspond to the observed outfiow velocities of Lyman 
break galaxies (LBGs) (e.g., Pettini et al. 2002). 

We have also implemented a prescription to insert the initial MHD radio jets associated 
with the SMBH formation, which in turn is based on the assumption that SMBH formation 
is directly related to merger events. Details will be presented elsewhere. 



3. Numerical Scheme for Ideal MHD 



3.1. Ideal MHD Equations 

Without expansion factor and external sources, the Equations ([I]) to (jl]) can be written 
in the conservative form as (in Cartesian coordinates): 



da OF dG OH 
where q = (p, pvi, pv2, pvs, Bi, B2, -B3, -E)*, and the fiux functions are 



G 
H 



( m 

pvj — Bf + p* 
PV1V2 - B1B2 
PV1V3 - B1B3 


-^2 



PV2 

PV2V1 - B2B1 

pvl - Bl + p* 

PV2V3 - B2B3 



-ill 



PV3 

PV3VI - B3BI 
PV3V2 - B3B2 

pvl -Bl + p* 

^2 





(7) 



\{E + p*)vi - Bi{B ■ {E + p*)v2 ~ B2{B ■ ^) {E + p*)v3 - BsiB ■ y) / 
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AU variables carry their usual meaning, 



is the total pressure, 

fll = V2B2, - V3B2, fl2 = V3B1 - V1B3, fls = V1B2 - V2B1 



(9) 



(10) 



are the electromotive force (EMF, defined via f2 = v x B). We also assume the ideal gas 
law for equation of state, and hence the total energy is 



^ = 2"^ +,-1 ' 2 



Prom the flux functions, we can obtain the Jacobian matrices 
A ^ , dF ^ , ^ dG ^ ^ ^ dH 



(12) 



Since G and H can be related to F through proper index permutation, Ay and Az 
should have similar structure as A^. The eigenvalues and eigenvectors for A^ have been 
extensively studied by many authors (see Brio & Wu 1988, Ryu & Jones 1995, etc.). Here 
we will write out the eigenvalues without explanation. In a direct extension of the ID system 
with magnetic field component Bi and total magnetic field B, the eigenvalues for Ax{q) are 



Xlj = Vi±Vf, X2,6 = Vi±Va, X3,5=Vi±Vs, X4 = Vi, 

where Va — \/B\l p is the Alfven speed based on B^-, and 



(13) 




B2 
P 



C2 + 



B^ 



— 4c^f ^ 



(14) 



(15) 



are the speeds of fast and slow magneto-sonic waves respectively, and Cg = Vtp/p the 
acoustic wave speed. 

In an alternative formulation, Powell et al. (1999) have used an 8 x 8 eigen-system. Its 
corresponding eigenvalues are 



Ai,8 = ±^;/, X2J^Vi±Va, X3^g^Vi±Vs, X4^5^Vi 



(16) 
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The corresponding eigenvectors have been given previously by many authors. We have 
adopted the eigenvector set given by Powell et al. (1999) in our code. 

For the dimensional split approach, the CFL condition along each direction can be 
calculated with (assume the same timestep is used) 



CFL = dt ■ min 



l^ll + Vf \V2\ + Vf \Vs\ + Vf 



dxi 



dxo 



dx^ 



(17) 



where the minimum is over all the cells. For the unsplit approach, the CFL condition is 
given as 



CFL ^ dt ■ min 



\Vi\ + Vf 

dxi 

where again, the minimum is over all the cells. 



+ 



\V2\ +'Vf ^ \V3\ + Vf 



dxo 



dxr^ 



(18) 



3.2. Numerical Solvers for Ideal MHD 

In recent years a variety of numerical algorithms for MHD based on the Godunov method 
have been developed. Early development of higher order Godunov schemes for MHD focused 
on interpreting the MHD equations as a simple system of conservation laws. This was 
done by Brio and Wu (1988), Zachary, Malagoh and Colella (1994), Powell (1994), Dai 
and Woodward (1994), Ryu and Jones (1995), Roe and Balsara (1996), etc. There are two 
important extensions to the basic HD algorithm for use of MHD. The first is an extension 
of the Riemann solver used to compute the fluxes of each conserved quantity to MHD; the 
second is to preserve the divergence-free constraints V • B = numerically. 

A typical second-order Godunov scheme for hyperbolic conservation law consists of three 
steps: 

1. Predict the time-centered cell- interface values, f/^"*"^/^ and U'^^^'^ . 

2. Using these Left and Right states at each cell-interface, solve the Riemann problem 
to determine the time-centered fluxes, F(C/^^^''^, C^^"'"^''^), possibly modified by the 
additional artificial viscosity. 

3. Compute the cell-centered conservative variable at the next time step using a conser- 
vative update with the time-centered fluxes. 

The first step is a predictor step via data reconstruction and time evolution. The second 
step is a Riemann solver. 
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3.2.1. Data Reconstruction 

We included several data reconstruction approaches in our code. We used piecewise 
linear reconstruction along with several slope limiters, including a general minmod, 

Aui = minmod (9{ui+i - Ui), ^(wj+i - Ui-i),9{ui - Ui-i) \ , (19) 



where Aui is a limited slope for cell [Xj^_i,x^_^i], and 9 E [1,2]. Note that a larger 9 
corresponds to less dissipation, but still non-oscillatory limiters. For ^ = 2, it becomes a 
Woodward limiter. These minmod-type limiters are not smooth functions with respect to u. 
Another limiter we used is van Albada-type, which is a smooth function of u, 

2{ui+i - Ui){ui - Ui^i) + e 1 



where e is a tiny positive constant in case of Ui = Wj+i and Ui — Ui^i. This limiter is smooth 
and it preserves the monotonicity of the original profile of u. 

For robustness, we have also implemented the characteristic limiting data reconstruc- 
tion, which applies the limited slope on wave-by-wave decomposition space obtained from 
local solutions of Riemann problems. In addition, we have also implemented the piecewise 
parabohc (Colella & Woodward 1984) and central WENO (Levy et al. 2000; 2002) re- 
construction schemes. For robustness, the data reconstruction is applied to the primitive 
variables rather than the conservative variables. 



3.2.2. Predictor for Time-centered Quantities 

We have implemented several approaches as the predictor for the time-centered quanti- 
ties. The simplest approach is the Hancock predictor, 

«r^ = < + ^^(F«+^A;<)-F«-iA:<)). 

where A^uf is the limited slope that is obtained from the data reconstruction. Then we use 
^ and AxUf to construct the left and right interface values at time ^""^2. In x direction, 
we have 

<^=<^^ + 2^' <i = ^^i' - 2^^- (21) 
Since the Hancock prediction does not need the characteristic decomposition, it can be done 
fast. Note that the limited slopes are reconstructed only once during the whole process. We 
do not need to calculate the hmited slope for intermediate state 
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We have also implemented the piecewise-linear method (PLM) via upwinding charac- 
teristic tracing (Colella and Glaz, 1985). 

3.2.3. Riemann solvers 

We implemented three different Riemann solvers, including Roe (Powell et al. 1999), 
HLL/HLLE [see Harten, Lax, & van Leer (1983); Einfeldt (1988); Einfeldt et al. (1991)], 
and HLLC (Li 2005). To improve the robustness, we also combined the HLLE and Roe's 
Riemann solves together in our code as proposed by Janhunen (2000). 

3.3. Dual Energy Formulation for High Mach and Low-/5 Plasma Flows 

In cosmological simulations, hypersonic ffows with a large Mach number (e.g., M > 100) 
appear to be very common, and they present a problem in our numerical schemes so far 
because the thermal energy is very small compared to the kinetic energy. The pressure, 
which is proportional to the thermal energy and used extensively in the numerical schemes, 
cannot be tracked accurately, because the discretization errors made in computing the total 
energy and the kinetic energy can be large enough to result in negative pressure. A similar 
problem arises in modeling the MHD flow with a low plasma beta {(3 = 2P/B^). 

Dynamically, this is not a large problem because the negative pressure can be replaced 
with a nominal floor value and the total energy budget of the flow remains unaffected. 
However, if the temperature distribution is considered for other reasons, (e.g., radiative 
processes), the thermal energy must be tracked accurately. 

As suggested by Ryu et al. (1993) and Bryan et al. (1995), we present dual energy 
formulations to track the thermal energy accurately. For comparison, we implement two 
approaches. The first is to solve the modified entropy Eq. ([5]), which is in conservative form 
and easy to implement for solvers that require eigen-decomposition. The other is to solve 
the internal energy density equation in Eq. ([6]). 

A similar set of Riemann solvers described in § fl3.2.3p are also implemented when solving 
the dual energy formulation. Since the internal energy equation is not in conservative form, 
we calculated only the flux V ■ (pev) in Eq. with the Riemann solver. The Roe's solver 
requires eigen-decomposition for a conservative system and cannot be applied to Eq. 
directly. For the modified entropy equation we adopt the eigen-decomposition from 
Balsara and Spicer (1999b) in our Roe's solver. The HLL and HLLC solvers can be applied 
to flux V ■ (pev) as well as to Eq. ([5]) in a straightforward manner. The pV ■ v in Eq. ([6]) 
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can be discretized as 



ra+1/2 n+1/2 
;+l/2 _ P^+l/2 +Pi-l/2 ( n+1/2 
~ ^ l^j+1/2 



— V 



n+1/2 
'i-1/2 



)/Aa; 



(22) 



where the time-centered values of j). 



n+1/2 
'i+1/2 



and V, 



n+1/2 
'j+1/2 



are obtained from the Riemann solver. 



It is necessary to use the total energy equation as much as possible for proper conversion 
of kinetic to thermal energy. When the total energy equation fails to give the positive 
pressure, we switch to solve the modified entropy or internal energy equation. This is done 
by monitoring the values of the internal energy e that comes from the total energy equation. 
When the ratio between the internal energy e and the total energy E satisfies e/E < rj, 
Eq. ([5]) or Eq. is solved and the total energy is updated with the new internal energy. 
We remark that as long as the parameter rj is small enough, the dual formulation will have 
no dynamical effect. After testing a large number of problems, we have chosen rj = 0.008. 
In Ryu et al. (1993), rj = 0.02 combined with a shock detection is used. In Bryan et al. 
(1995), T] = 0.001 combined with a shock detection is used. In our code, we tested t] = 0.008 
with and without shock detection. Although both Eqs. ([5]) and will not give a correct 
solution for the shock, we have not found any problem fails using e/E < t] only without 
shock detection. On the other hand, if e/i? < ?7 is used with a shock detection scheme, the 
internal energy of the cell near the shock usually is set to a nominal floor value, which is 
incorrect. 



Maintaining the divergence-free condition V ■ B = is important in MHD simulations. 
Brackbill and Barnes (1980) find that the non-zero divergence, caused by discretization 
errors, can grow exponentially during the computation and destroy the correctness of the 
solutions. Hence they proposed a divergence-cleaning approach, which solves an extra global 
elliptic (Poisson) equation to recover V ■ B = at each time-step. Balsara & Kim (2004) 
find that the divergence-cleaning is significantly inadequate when used for many astrophysical 
applications. The constrained transport (CT) approach uses a staggered grid, which places 
the magnetic field variables at the cell-face to keep V • B to the accuracy of machine round- 
off error. Recently, this approach has been combined with various modern shock-capturing 
algorithms by many authors (see Toth 2000 and its references). Toth (2000) proposed a 
constrained transport/central-difference (CT/CD) method that uses cell-centered values of 
the magnetic fields. This method works well with single grid but it is unknown currently 
how to apply it to adaptive grids. In the current stage of CosmoMHD, we deal with only 
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a single grid. Therefore, we adopted the flux-CT method of Toth (2000) to preserve the 
divergence-free condition. 



4. Numerical Tests 

4.1. Shock-tube Test 

We first tested the code with an MHD shock tube problem on a two-dimensional (2D) 
Cartesian grid. The shock tube is usually posed as a ID problem, and it is a standard test 
for a numerical scheme to handle the discontinuities, e.g., shock and contact. Solving it on a 
2D grid can verify if a multi-dimensional scheme works properly and recovers the ID solution 
profile. 

We chose an example that has been used by Toth (2000) to compare several numerical 
schemes for MHD. We adopted the same initial and boundary conditions as Toth (2000). 
The initial left and right states are 



rn„ „ r, R R 1(1.08,1.2, 0.01, 0.5, 0.95, 2/v^,3.6/v^,2/v^), left, 

(p, , v^, p, i^ii, i^^, B,) - I 4/v/4^, 2/ V4^), right. 



where || refers to the direction along the normal of the shock front, ± refers to the direction 
perpendicular to the normal of the shock front but still in the computational plane, and 
z refers to the direction out of the plane. Note that Vz and Bz components are non-zero. 
Therefore, it is often called the 2.5D MHD shock tube problem. We first solve this problem 
using a ID grid and then solve it as a 2D problem with a rotation angle a = 45° between 
the shock interface and y-axes. The same boundary conditions as Toth (2000) and 256 cells 
in each direction are used. The computation is stopped at t = 0.1/ \/2 before the fast shocks 
reach the left and right boundaries. The numerical solution is compared with a reference 
solution that is calculated with a ID grid and 1600 cells. 

For this example, only the total energy equation is solved. The modified entropy equa- 
tion and internal energy equation are not used. Fig. [1] shows the results at the final time. 
The parallel component of the magnetic field is preserved very well by our CT scheme, which 
can preserve the divergence-free condition (V • B = 0) to machine round-off error. 



4.2. One-dimensional MHD Caustics 



This example is taken from Ryu et al. (1993). It was proposed as a test for hydro 
cosmology code to handle the thermal energy correctly. The formation of ID caustics has 
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Fig. 1. — Plots of density, magnetic components in parallel and vertical direction along the 
normal of the shock front, and pressure of the 2.5D shock-tube problem. The reference 
solution is calculated with ID grid and 1600 grid cells. The ID and 2D results are calculated 
with 256 and 256x256 cells respectively. Output is at t=0.1 for the reference solution and 
ID result, and at t = 0.1/^/2 for 2D results with shock angle of « = 45°. The parallel 
component of the magnetic fields {B\\) is conserved with an accuracy of machine round-off 
errors. 
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been calculated from the initial sinusoidal velocity field along the x-direction with the peak 
value l/27r. The initial density and pressure has been set to be uniform with p = 1 and 
p = 10"^'^. Hence the initial peak velocity corresponds to a Mach number of 1.2x10^. In the 
simulations, the pressure can easily become negative for this type of high Mach flow. 

We test this problem with and without magnetic fields. If the magnetic fields are present, 
we set = Bz = 0, and only By is nonzero. Different values of By are tested. Fig. [2] shows 
the results at t = 3 with 1024 cells. For this problem, the total energy version of the code 
will give a wrong pressure floor value, which is not shown here. The results shown in Fig. 
[2] are calculated using the combination of the total energy version with a modifled entropy 
equation. It is clear that the strong magnetic flelds (5=0.05) can suppress the density peak, 
while relatively weak magnetic flelds (i?=lE-3) have little impact on the density and pressure 
proflle. 
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Fig. 2. — Plots of density, gas pressure, total pressure, and magnetic component {By) of ID 
MHD caustics test at t = 3. The modifled entropy equation is solved. 
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We also solve the internal energy equation to follow the adiabatic changes correctly in 
the pre-shock region. Figure [3] shows the comparison between the modified entropy version 
and the internal energy version for By = 0.05 case. The results are almost identical. Figure 
E] also shows the comparison with different Riemann solvers. Apparently the HLL solver is 
very diffusive in the post-shock region so that it cannot recover the gas pressure correctly in 
the middle region. 

Fig. 3. — Comparison of different approaches (left) and Riemann solvers (right) with the 
plots of gas pressure of ID MHD caustics test at t = 3. Initial By = 0.05. 



4.3. One-dimensional Cosmological Pancake Collapse 

This example is also taken from Ryu et al. (1993). We set up the problem as a ID 
pancake (i.e., the Zel'dovich pancake problem) in a purely baryonic Universe with Q = 1 and 
h = 0.5 in the comoving coordinates. The MHD solver with self-gravity module is tested. 
Initially, at = 1, which corresponds to Zi = 20 in this test, a sinusoidal velocity field with 
the peak value 0.65/(1 + Zi) in the normalized units has been imposed in a box with the 
comoving size 64/i~^Mpc. The baryonic density and pressure have been set to be uniform 
with {p,p) = (1,6.2 X 10^^) in the normalized units. The calculation has been done with a 
different number of cells, with and without magnetic fields. When the magnetic fields are 
present, we use only the By component, and set B^ = B^ = 0. 

Figure m shows the results with By = 0, i.e., without the magnetic fields. The results for 
both the modified entropy version and internal energy version are shown. For comparison, we 
also show the results of TVD-ES code (Ryu et al. (1993)), running with the mass-diffusive 
correction. 

Figure E] shows the results with magnetic fields. Different magnetic fields have different 
impact on the shock location and the peaks of density and magnetic fields. With a small 
By = 1.3e — 6, the magnetic fields have little impact on the shock and the density profile. 
The magnetic fields collapse in the same way as the density field does. As we increase 
the magnitude of the magnetic fields, the density peak in the post-shock region and shock 
location show large changes, while the profile in the pre-shock region has little change. 

Figure [6] shows the comparison for different Riemann solvers applied to the modified 
entropy equation. The results for HLLC and the more sophisticated Roe's solver are almost 
identical, while the result for HLL solver is not so good at the post-shock region. 
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Fig. 4. — One-dimensional gravitational collapse of a pancake in the comoving coordinates 
with different codes. The number of cells is 1024. No magnetic fields are apphed {B — 0). 
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Fig. 5. — One-dimensional gravitational collapse of a pancake in the comoving coordinates 
with different magnitude of magnetic fields. The number of cells is 1024. 
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Fig. 6. — One-dimensional gravitational collapse of a pancake in the comoving coordinates 
with different Riemann solvers. The number of cells is 1024. By = 10~^. 

4.4. Comparison with Other Codes for 3D Adiabatic CDM Universe 

We have carried out a comparison of our cosmological MHD code with the original 
TVD-ES code of Ryu et al. (1993). We computed the adiabatic evolution of a purely 
baryonic Universe but with an initial CDM power spectrum with the following parameters: 
Q = Qh = 1, h = 0.5,, n = 1 and (Xg = 1, and the size of the computational cube is 
L = 64/i^^Mpc. We use the Bardeen et al. (1986) transfer function to calculate the power 
spectrum of the initial density fluctuations. This test problem is identical to that of Kang 
et al. (1994) except that the random numbers used to generate the initial condition are 
different. The Universe is evolved from z = 30 to z = 0. We use 256^ cells for each 
simulation performed. The comparisons are made at the final epoch, z = 0. 

Figure [7] shows a comparison of a mass- weighted histogram of the baryonic temperature 
at z = without magnetic fields between the original TVD-ES code and CosmoMHD. The 
mass-weighted histogram is calculated by using the temperature as x-axis and total mass in 
each temperature bin as the y-axis. They are in a good agreeement. 

We have also performed another run with initial magnetic fields, 3^ = 3^ = 0, By=2.'oE- 
9 Gauss (which corresponds to 4.3176E-7 in code unit). Figure [H] shows a shce of the density 
and magnetic energy at 2 = 0. We see that the distributions of density and magnetic field 
are strongly correlated, as expected. To demonstrate the impact of the CT method, we 
set the initial i?j,=lE-5 Gauss, which corresponds to 0.0017 in code unit. Figure [9] shows 
the comparison of the divergence of the magnetic fields, averaged over the entire box, as a 
function of redshift. It is clear that without the CT method, the divergence of the magnetic 
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fields grows in tandem with the growth of the magnetic field strength, and simulation results 
would hardly be meaningful. Since single-precision is used in our whole simulation, the 
divergence from the CT method is close to the round-off error. The temperature plots 
are shown in Fig. (TUl The plasma beta becomes very small (10~^) in many regions with 
these large magnetic fields. The dual-formulation must be used to track the internal energy 
accurately. 




Fig. 7. — Mass-weighted temperature histogram at present (redshift z=0) for the 3D simu- 
lations of a purely baryonic adiabatic Universe. 



Fig. 8. — Density and magnetic energy at present for the 3D simulations of a purely baryonic 
adiabatic Universe. Left: log density plot at z=0.5 slice. Right: log magnetic energy plot at 
z=0.5 slice. Initial magnetic fields are = = and By=A.3176E-7. 
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Fig. 9. — The scaled divergence (/iV-B/5) of the magnetic fields with and without constraint 
transport (CT) for the three-dimensional simulations of a purely baryonic adiabatic Universe 
during the whole time history. The initial By at Z = 30 is 0.0017. 



Fig. 10. — Temperature plots of a slice at present for the 3D simulations of a purely baryonic 
adiabatic Universe. Left: log temperature plot without using CT at z=0.5 slice. Right: log 
temperature plot with CT at z=0.5 slice. Initial magnetic fields are = Bz = Q and 



5^^=0.0017. 
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5. Discussion and Conclusions 

We have developed and tested a modern cosmological MHD code (CosmoMHD). In 
the CosmoMHD code, five sets of equations are evolved simultaneously: (1) the ideal MHD 
equations for magneto-gasdynamics, (2) rate equations for multiple species of different ioniza- 
tional states (including hydrogen, helium and oxygen), (3) the Vlasov equation for dynamics 
of coUisionless particles (including dark matter particles and stellar particles), (4) the Pois- 
son's equation for obtaining the gravitational potential field and (5) the equation governing 
the evolution of the intergalactic ionizing radiation field, all in cosmological comoving coor- 
dinates. The MHD solver consists of several high-resolution schemes with finite-volume and 
finite-difference methods. The divergence-free condition of the magnetic field is preserved to 
an accuracy due just to round-off errors. Additional implemented physical processes include: 

(1) detailed cooling and heating processes due to all the principal line and continuum pro- 
cesses for a plasma of primordial composition as well as coohng and heating due to metals, 

(2) star formation processes and feedback processes from star formation including UV radi- 
ation and galactic superwinds, and (3) formation of supermassive black holes and associated 
feedback processes including radio jets/lobes. 

The CosmoMHD code developed here may be applied to a variety of astrophysical/cosmological 
problems with or without self-gravity. We intend to study the formation and evolution of 
galaxies, clusters of galaxies and the intergalactic medium using this accurate code. Most 
importantly, we would like to investigate the role played by the AGN feedback in the form of 
powerful radio jets/lobes in regulating structure formation on scales from the cores of clus- 
ters of galaxies (~ lOkpc) through giant lobes (~ lOOkpc— IMpc) to large scales (> IMpc). 
It seems prudent to check if many fundamental cosmological applications, such as weak grav- 
itational lensing, formation and evolution of X-ray/SZ clusters and matter power spectrum 
with baryonic oscillation features, may be affected by this feedback effect. A thorough un- 
derstanding of AGN feedback may be necessary to reduce systematic errors to levels that are 
commensurate with statistical errors that many of these important observations are expected 
to be able to achieve. 

We thank Eric Johnson for coding a prescription to form initial realistic MHD radio 
jets associated with SMBH formation. This work was carried out under the auspices of the 
National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos 
National Laboratory under Contract No. DE-AC52-06NA25396, and is supported by the 
Laboratory Directed Research and Development programs at LANL. It is also supported in 
part by grants AST-0407176, AST-0507521, NNG05GK10G and NNG06GI09G. 
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